# 加载相关包
library(openxlsx)
library(geepack)

# 读取数据
data_test = read.xlsx("./output_data1.xlsx", 3)
data_anal = data_test[, c(1, 9:13, 14:25, 27, 29:33)]
varsx = colnames(data_anal)[-c(1:6, 29:33)]
varsy = colnames(data_anal)[c(2:6)]
varsc = colnames(data_anal)[c(29:33)]

# 转换变量类型：因变量为数值、自变量为因子
data_anal[varsy] = lapply(data_anal[varsy], as.numeric)
data_anal[varsx] = lapply(data_anal[varsx], as.factor)

wb = createWorkbook();
# for循环将分析结果输出
# 外循环控制因变量的个数
for (j in 1:length(varsy)) {
  out = NULL
  extract_data = data_anal[complete.cases(data_anal[, j + 1]), ]           # 去掉因变量的空行
  # 内循环控制自变量的个数
  for (i in 1:length(varsx)) {
    extract_data = extract_data[complete.cases(extract_data[, i + 6]), ]   # 去掉自变量的空行
    fit = geeglm(formula = as.formula(paste(varsy[j], "~", varsx[i])),     # 使用 as.formula() 生成公式
                 data = extract_data,
                 id = id, family = "binomial",
                 corstr = "ar1", scale.fix = TRUE)
    out1 = summary(fit)$coefficients
    out1$variable = varsx[i]
    # 追加结果：内循环第一次执行时 out 为空，会直接赋值为 out1；以后都是追加行
    out = if (is.null(out)) out1 else rbind(out, out1)
  }
  # 输出为单个表格
  addWorksheet(wb, sheet = varsy[j])
  writeDataTable(wb, sheet = varsy[j], x = out, rowNames = T)
}

saveWorkbook(wb, "./output.xlsx", overwrite = T)